%----------------------------------------------
% Formatting and file path
%----------------------------------------------

clear;
clc;
format bank; % Set formatting to display numbers without e format

%----------------------------------------------
% Read in data
%----------------------------------------------
Folder = cd; 
Folder = fullfile(Folder, '../../Data/structural/');
cd(Folder); 
for j = 3:5 
    %input = append('intermediate/data_calib_w', string(j), '.csv');
    input = ['intermediate/data_calib_w', num2str(j), '.csv'];
    disp(input)
    data = csvread(input);
    hhid=data(:,1);
    l_prev_w=data(:,2);
    h_prev_w=data(:,3);
    tot_hrs_prev_w=data(:,4);
    k_prev_w_post_trans=data(:,5);
    hired_l_prev_w=data(:,6);
    w_prev_w=data(:,7);
    case_prev_w=data(:,8);
    l_wj=data(:,9);
    h_wj=data(:,10);
    tot_hrs_wj=data(:,11);
    k_wj=data(:,12);
    hired_l_wj=data(:,13);
    w_wj=data(:,14);
    case_wj=data(:,15);
    k_prev_w_pre_trans=data(:,16);
    H=data(:,17);

    N=length(hhid);    % Number of HHs
    R=3650;            % Time endowment constraint
    Hire_in_max=1400;  % Maximum number of hours that can be hired in

    a=-1.88e-10;       % Estimated parameters
    b=0.0009368;
    beta=0.3499541;

    uns_threshold=10380;

    phi = 0.7;          % Ratio of hired in labour to market wage
    rho = 0;            % Share of capital profit from liquidating


    %----------------------------------------------
    % Simplify terms
    %----------------------------------------------

    fk_wj = a.*k_wj.^2 + b.*k_wj;
    fk_prev_w = a.*k_prev_w_pre_trans.^2 + b.*k_prev_w_pre_trans;

    %----------------------------------------------
    % Calibrate psi and A parameters
    %----------------------------------------------

    psi_h = zeros(N,1);
    psi_l = zeros(N,1);
    A = zeros(N,1);

    for i = 1:N
        % Mixing with hired-in labour in Year 2
        if case_wj(i) == 1
            A(i) = phi * w_wj(i) / (fk_wj(i)*beta*(l_wj(i) + hired_l_wj(i))^(beta-1));
            psi_h(i) = w_wj(i) / (phi*l_wj(i) + h_wj(i));
            psi_l(i) = phi^2 * psi_h(i);
        end

        % Livestock rearing only with hired-in labour in Year 2
        if case_wj(i) == 3
            if case_prev_w(i) == 6
                psi_h(i) = w_prev_w(i)/h_prev_w(i);
                A(i) = phi * w_wj(i) / (fk_wj(i)*beta*(l_wj(i) + hired_l_wj(i))^(beta-1));
                psi_l(i) = phi * w_wj(i) / l_wj(i);
            elseif case_prev_w(i) == 7
                A(i) = phi * w_wj(i) / (fk_wj(i)*beta*(l_wj(i) + hired_l_wj(i))^(beta-1));
                psi_l(i) = phi * w_wj(i) / l_wj(i);
            end
        end

        % Mixing without hired-in labour in Year 2
        if case_wj(i) == 2
            if case_prev_w(i) == 6
                psi_h(i) = w_prev_w(i)/h_prev_w(i);
                psi_l(i) = (w_wj(i) - psi_h(i)*h_wj(i))^2/(psi_h(i)*l_wj(i)^2);
                A(i) = (psi_l(i)*l_wj(i) + sqrt(psi_h(i)*psi_l(i))*h_wj(i))/(fk_wj(i)*beta*l_wj(i)^(beta-1));
            end
        end
    end

    % Replace missing psi_h with maximum

    psi_h(psi_h==0) = max(psi_h);

    % Rest of mixing without hired-in labour in Year 2
    for i = 1:N
        if case_wj(i) == 2
            if case_prev_w(i) == 7
                psi_l(i) = (w_wj(i) - psi_h(i)*h_wj(i))^2/(psi_h(i)*l_wj(i)^2);
                A(i) = (psi_l(i)*l_wj(i) + sqrt(psi_h(i)*psi_l(i))*h_wj(i))/(fk_wj(i)*beta*l_wj(i)^(beta-1));
            end
        end
    end
    
    T = table(hhid, A, psi_h, psi_l);
    %filename=append('output/w', string(j),'parameters.xlsx');
    filename = ['output/w', num2str(j), 'parameters.xlsx'];
    writetable(T,filename);
end